{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Thermo Properties\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reading Polynomial parameters and Calculating Thermodynamic parameters\n",
    "\n",
    "The thermo class given will be used to calculate the thermodynamic parameters. Also, you'll need [thermoDataNASA-9.yaml].  This file contains two arrays for each species for two separate temperature ranges: a_lo is for $200<T<1000$ K and a_hi is for $1000<T<6000$ K. We have followed Cantera nomenclature (http://www.cantera.org/docs/sphinx/html/cti/species.html#the-nasa-9-coefficient-polynomial-parameterization).\n",
    "\n",
    "\n",
    "$$\\frac{C^0_p(T)}{R_u}=a_0T^{−2}+a_1T^{−1}+a_2+a_3T+a_4T^2+a_5T^3+a_6T^4$$\n",
    "\n",
    "$$\\frac{H^0(T)}{R_u T}=−a_0T^{−2}+a_1 \\frac{lnT}{T}+a_2+\\frac{a_3}{2}T+\\frac{a_4}{3}T^2+\\frac{a_5}{4}T^3+\\frac{a_6}{5}T^4+\\frac{a_7}{T}$$\n",
    "\n",
    "$$\\frac{s^0(T)}{R_u}=−\\frac{a_0}{2}T^{−2}−a_1T^{−1}+a_2 lnT+a_3T+\\frac{a_4}{2}T^2+\\frac{a_5}{3}T^3+\\frac{a_6}{4}T^4+a_8$$\n",
    "\n",
    "where $R_u = 8.31446 J/mol/K$ is the universal gas constant.\n",
    "\n",
    "Note that Gordon and McBride start polynomial coefficient numbering at 1 and the last coeffients for enthalpy and entropy are named respectively $a_7 = b_1$ and $a_8 = b_2$. (https://ntrs.nasa.gov/api/citations/20020085330/downloads/20020085330.pdf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This cell of code give you an implementation of the above formulas.\n",
    "import yaml\n",
    "import numpy as np\n",
    "\n",
    "class thermo:\n",
    "    \n",
    "    def __init__(self, species, MW) :\n",
    "        \"\"\"\n",
    "        species: input string name of species in thermoData.yaml\n",
    "        M: input (species molecular weight, kg/kmol)\n",
    "        \"\"\"\n",
    "        \n",
    "        self.Rgas = 8.31446      # J/mol*K\n",
    "        self.M    = MW\n",
    "    \n",
    "        with open(\"thermoDataNASA-9.yaml\") as yfile:\n",
    "           yfile = yaml.safe_load(yfile)\n",
    "        \n",
    "        self.a_lo = yfile[species][\"a_lo\"]\n",
    "        self.a_hi = yfile[species][\"a_hi\"]\n",
    "        \n",
    "        self.T_lo = 200.\n",
    "        self.T_mid = 1000.\n",
    "        self.T_hi = 6000.\n",
    "        \n",
    "    def cymalalo(self) :\n",
    "        return self.a_lo\n",
    "\n",
    "    def cymalahi(self) :\n",
    "        return self.a_hi\n",
    "\n",
    "    def cp_mole(self,T) :\n",
    "        \"\"\"\n",
    "        return calorific value at cst p in units J/mol/K\n",
    "        T: input (K)\n",
    "        \"\"\"\n",
    "        if T<=self.T_mid and T>=self.T_lo :\n",
    "            a = self.a_lo\n",
    "        elif T>self.T_mid and T<=self.T_hi :\n",
    "            a = self.a_hi\n",
    "        else :\n",
    "            print (\"ERROR: temperature is out of range\")\n",
    "\n",
    "        cp = a[0]/T**2 + a[1]/T + a[2] + a[3]*T + a[4]*T**2.0 + a[5]*T**3 + a[6]*T**4\n",
    "        \n",
    "        return cp * self.Rgas\n",
    "        \n",
    "    #--------------------------------------------------------\n",
    "\n",
    "    def cp_mass(self,T) :\n",
    "        \"\"\"\n",
    "        return calorific value at cst P in units of kJ/kg/K\n",
    "        T: input (K)\n",
    "        \"\"\"\n",
    "        return self.cp_mole(T)/self.M\n",
    "\n",
    "    #--------------------------------------------------------\n",
    "    \n",
    "    def h_mole(self,T) :\n",
    "        \"\"\"\n",
    "        return enthalpy in units of J/mol\n",
    "        T: input (K)\n",
    "        \"\"\"\n",
    "        if T<=self.T_mid and T>=self.T_lo :\n",
    "            a = self.a_lo\n",
    "        elif T>self.T_mid and T<=self.T_hi :\n",
    "            a = self.a_hi\n",
    "        else :\n",
    "            print (\"ERROR: temperature is out of range\")\n",
    "\n",
    "        hrt = -a[0]/T**2 + a[1]*np.log(T)/T + a[2] + a[3]/2*T + a[4]/3*T**2.0 + a[5]/4*T**3 + a[6]/5*T**4 + a[7]/T\n",
    "        \n",
    "        return hrt * self.Rgas * T\n",
    "        \n",
    "    #--------------------------------------------------------\n",
    "\n",
    "    def h_mass(self,T) :\n",
    "        \"\"\"\n",
    "        return enthalpy in units of kJ/kg\n",
    "        T: input (K)\n",
    "        \"\"\"\n",
    "        return self.h_mole(T)/self.M\n",
    "\n",
    "    #--------------------------------------------------------\n",
    "        \n",
    "    def s_mole(self,T) :\n",
    "        \"\"\"\n",
    "        return entropy in units of J/mol/K\n",
    "        T: input (K)\n",
    "        \"\"\"\n",
    "        if T<=self.T_mid and T>=self.T_lo :\n",
    "            a = self.a_lo\n",
    "        elif T>self.T_mid and T<=self.T_hi :\n",
    "            a = self.a_hi\n",
    "        else :\n",
    "            print (\"ERROR: temperature is out of range\")\n",
    "        \n",
    "        sr = -a[0]/2/T**2 - a[1]/T + a[2]*np.log(T) + a[3]*T + a[4]/2.0*T**2.0 + a[5]/3.0*T**3.0+ a[6]/4.0*T**4.0+ a[8]\n",
    "        \n",
    "        return sr * self.Rgas\n",
    "        \n",
    "    #--------------------------------------------------------\n",
    "\n",
    "    def s_mass(self,T) :\n",
    "        \"\"\"\n",
    "        return entropy in units of kJ/kg/K\n",
    "        T: input (K)\n",
    "        \"\"\"\n",
    "        return self.s_mole(T)/self.M"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Entropy of  O2  at  298 K:  205.13344991630754  J/mol\n",
      "Enthalpy of  O2  at  298 K:  -4.406684233742955  J/mol\n",
      "Specific heat of  O2  at  298 K:  29.377442521060257  J/mol\n"
     ]
    }
   ],
   "source": [
    "# Example on how to use the class thermo to calculate thermodynamics \n",
    "# quantities for a specific species.\n",
    "Molecule = \"O2\"    \n",
    "MW_molecule = 32               # kg/kmol\n",
    "\n",
    "t = thermo(Molecule,MW_molecule) # thermo object;\n",
    "Temperature = 298       #K\n",
    "print('Entropy of ',Molecule,' at ',str(Temperature), 'K: ',t.s_mole(Temperature),' J/mol')\n",
    "print('Enthalpy of ',Molecule,' at ',str(Temperature), 'K: ',t.h_mole(Temperature),' J/mol')\n",
    "print('Specific heat of ',Molecule,' at ',str(Temperature), 'K: ',t.cp_mole(Temperature),' J/mol')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# If you wish to use this piece of code embedded in another code, you might need the following lines:\n",
    "# Careful, you might need the following function the first time you are using # the code: \"!pip install ipynb\"\n",
    "from ipynb.fs.full.ThermoPropertiesNASA import thermo\n",
    "\n",
    "# Try it out in a new script, if it works, it should display the same results than the cell above\n",
    "# If it works, just copy the content of the previous cell in your new script, and adapt it to your needs\n"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
